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Investigating ultra-long gravitational waves with 
measurements of pulsars rotational parameters 
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ABSTRACT 

A method is suggested to explore the gravitational wave background (GWB) in the 
frequency range from 1CP 12 to 10~ 8 Hz. That method is based on the precise mea- 
surements of pulsars' rotational parameters: the influence of the gravitational waves 
(GW) in the range will affect them and therefore some conclusions about energy den- 
sity of the GWB can be made using analysis of the derivatives of pulsars' rotational 
frequency. The calculated values of the second derivative from a number of pulsars 
limit the density of GWB Q, gw as follows: tt gw < 2 x 10~ 6 . Also, the time series of 
the frequency v of different pulsars in pulsar array can be cross-correlated pairwise in 
the same manner as in anomalous residuals analysis thus providing the possibility of 
GWB detection in ultra-low frequency range. 

Key words: gravitational waves-pulsars: general-cosmology: miscellaneous- 
methods: data analysis 



1 INTRODUCTION 

Search for gravitational waves is one of the most important 
tasks in modern astronomy and physics. Various techniques 
are used to look for gravita tional waves in a very broad range 
of fre quencies fro m 10~ 18 ( Baskara n. Grishchuk fc Polnarevl 
120061 ') to 10 3 Hz l|LIGOll2008D ~ 

Pulsar timing provides a unique access for obser vations 
in a low-frequency band (10~ 7 Hz < f aw < 10~ 9 H z) (|Sazhinl 
ll978l : lDetweiienll979l : lB"ertotti. Carr fc Reeslll983l ). 

Propagation of pulsar signal in space-time perturbed by 
a stochastic gravitational wave field results in apparent de- 
viations of pulsar rotational frequency; that influence can 
be sought in anomalous residuals of pulses arrival times. 
RMS of that residuals can be transferred to the upper lim- 
its on density of gravitational wave background (GWB). 
That background affects signals from all pulsars, and com- 
mon response c an be extracted from c e rtain correlations 
in timing series dHellings fc Downs! 1 19831 ; IJenet et al. I l2005l : 
lAnholm et al.l l2008). thus pulsar timing can not only place 
upper limits, but can also detect the presence of gravita- 
tional waves as well. Like any other gravitational wave de- 
tector, pulsar timing has its own frequency limitations. The 
lowest frequency of the gravitational wave that can be ob- 
served corresponds to T~^, where T a ^ s - total time span 
of pulsar observations (usually years). Influence of gravita- 
tional waves of lower frequencies simply redefines observed 
values of derivatives of pulsar rotational parameters v, v. 
Sensitivity of the method decreases with increase of GW 
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frequency, so the highest attainable frequency is of order of 
10 -7 Hz. Less stringent limitations in the frequency region 
~ (cD) -1 < f gw < TT , where D - distance to the pulsar, 
can be put using pulsar in binary system as a precise clock. 
In that case, deviations from values of pulsar orbital param- 
eters predicted by General Relativity are treated as manifes- 
tations of GWB induced effect and th erefore some upper lim- 
its on GWB density can be obtained |Bertotti. Carr fc Reesl 
1 19831 : lKopeikinlll997l ; iPotapov et al.ll2003l ). That method ex- 
tends achievable frequency region down to 10~ 12 Hz. In this 
paper a similar method is used: gravitational waves affect 
the derivatives of pulsars rotational frequency - these values 
absorb the effect caused by gravitational waves of ultra-low 
frequencies (~ (cD) _1 < f gw < T~^ s ). But the information 
about gravitational wave background is not lost completely 
and analysis of rotational parameters can constrain the char- 
acteristic s of GWB in the frequency range from 10 -12 to 
10" 8 Hz (|Bertotti. Carr fc Reesl Il983l ). That method com- 
plements the one that is based on observations of pulsar 
binaries. Moreover, using time series of v(t) in the same 
manner as time series of ToA residuals (IJenet et al.ll2005l ) 
and looking for the correlations between them we can try to 
detect GWB in the ultra- low frequency region. 

Here we briefly review possible sources of GW in that 
frequency range. The gravitational radiation from astro- 
physi cal sources like SMBH -binaries are negligible in the re- 
gion (|Jaffe fc Backerll2003l ). so the main aim for the searches 
is GWB of cosmological origin. There are a lot of possi- 
ble candidates for GW sources: relic gravitational waves 
from primordial fluctuations that were amp lified during the 
inflationary stage in the Early Universe l|Maggiorel l200Cl : 
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iGrishchuk et al.ll200ll ). gravita tional waves f rom phase tran- 
sitions in early Universe (e.g. dWittenl 1984)). gravitationa l 
waves from the string network ( Damour fc Vilenkinl f2005). 
The spectrum of the GWB depends on its sources, it can be 
either flat (Harrison-Zeldovich spectrum) or it can possess 
different features (e.g. peak a t ~ 10~ 12 Hz in some mod els 
of GWB produced by strings (|Vilenkin fc Shellard[l994l )). 



2 CONSTRAINTS FROM THE SECOND 
DERIVATIVE OF THE ROTATION 
FREQUENCY 

Precise pulsar timing allows us to measure seco nd derivative 
of rot ational frequency for a number of pulsars (|Hobbs et al.l 
120041 '). However, that value has a little in common with phys- 
ical value v that somebody would expect from any adopted 
laws of pulsar spin-down and usually surpasses the latter by 
several orders of magnitude. So, it is very plausible to regard 
the calculated value of v as caused by unknown factors that 
are intrinsic or extrinsic to pulsar. 

The technique o f calculation in the paper follows 
l|Baskaran et al.ll2008T l with parameter e = 0. We work in 
the framework of a slightly perturbed Minkowski space-time 
with coordinates x^ = (ct, x') and the metric given by: 



We can write down the final result: 



ds 



-c dt + (Sij + hij) dx 1 dx j 



(1) 



where hij is the gravitational wave perturbation. Firstly, 
we consider the simplest case of a single monochromatic 
wave and then generalize it to the case of a stochastic 
GWB. For a monochromatic grav itational wave the metric 
perturbation hg takes the form (jLandau fc Lifshit3 1 19751 : 
IMisner. Thorne fc Wheeled Il973h : 



hg = h Pije tk " x = h pge 



— i^kQct — k^x^^ 



(2) 



where h is the amplitude of the gravitational wave, k^ = 
(fco, ki) is the wave vector, and pik is the polarization tensor 
of the gravitational wave. We can introduce a set of two 
mutually orthogonal unit vectors h and mi orthogonal to 
the wave vector k{\ the polarization tensor pik has the form 
jLandau fc Lifshitzlll975l : IMisner. Thorne fc Wheeled] 19731) . 
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where vo is the unperturbed pulsar frequency in the absence 
of gravitational waves and 8v(t) = v(t) — vo is the variation 
of pulsar frequency due to their presence. D is the distance 
from the pulsar to the observer, e 1 is the unit vector tangent 
along this path (i.e. unit vector in the direction from pulsar 
to the observer). 

Gravitational wave will also affect derivatives of pulsar 
frequency analogously: 
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Our analysis can be generalized to a stochastic grav- 
itational wave background case. That background can be 
decomposed into spatial Fourier modes: 



hij(t,x l ) = 



l ) = f d :i k [hs{k\t)Pij (k')e ikixi +, 

** 1 n 



(8) 



where d k denotes the integration over all possible wave vec- 
tors, "c.c" stands for "complex conjugate" and s = 1,2 cor- 
responds to two linearly independent modes of polarization 
satisfying the orthogonality condition 



PiaP 



(9) 



The mode function h s (k l , t) corresponds to plane monochro- 
matic waves 



h s (k l ,t) = h s {k l ) e 



(10) 



Due to the linearity of the equations, shifts in derivatives of 
frequency can be presented as integrals: 



= Jd 3 kJ2 [hs(k')R +1 , + 2(t;k\s)+c.c], (11) 



Pik — -(U ± mi) (h ± rrik) 



(3) 



where ± corresponds to two independent states of circular 
polarization. Due to the transverse and traceless nature of 
gravitational waves, the polarization tensor satisfies the fol- 
lowing conditions 



Pikk' = 0, PikS'' 



0. 



(4) 



ft is convenient to introduce the wavenumber k = 
(Sijtfki) ' ,and a unit vector in the direction of wave prop- 
agation fc* = k % /k. The wavelength of the gravitational wave 
is related to the wavenumber by the equality k — 2n/X gm 
and the frequency of the gravitational wave f gw is related to 
the time component of the wave vector through the relation 
fco = 2irf aw /c. 

The effect of a gravitational wave upon the measured 
frequency of pu lsar signal was considered in (|Sazhinlll978l : 
|Petweiler|[l979l ). 



where indices +1, +2 refer to the first and the second deriva- 
tive respectively. Applying the results from the consideration 
of a single monochromatic wave, we can write down the con- 
tribution from a single Fourier component R+^^it; k 1 , s): 



R+i(t; k\ s) = 



R+2(t; k\ s) = 



-ike 
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(13) 



where the tilde over .R+1,+2 in the above expressions is 
introduced to indicate explicit factoring out of the gravita- 
tional wave amplitude ft. 

Usually we know only statistical properties of the grav- 
itational wave field. Stationary statistically homogeneous 
and isotropic gravitational wave field possesses the following 
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properties: 



<h s (k l ) h* 3 r(k' i ) >^^\s sa ,S 3 (k 



< h s (k l ) >= 0, 



(14) 
(15) 



167rfc 3 

where the brackets denote ensemble averaging over all possi- 
ble realizations, and Ph (k) is the metric power spectrum per 
logarithmic interval of k. Using (|15[l . we can calculate the 
statistical properties of the corresponding shifts in frequency 
derivatives i>{t) and v(t). Using pip and (|15p . and taking 
into account the orthogonality property ([9|, after straight 
forward calculations, we arrive at the following statistical 
expressions: 



5u(t) 
< — — > 
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dk 

^P h (k)R 2 +1 (k), 



dk 

^P h (k)R 2 +2 (k), 



£+i,a(fc) 



where we have introduced the transfer functions 



fl + i,a(t; k',s) 



(16) 
(17) 
(18) 
(19) 

(20) 



In the above expression dfl represents integration over 
the possible directions of gravitational wave (i.e. d k = 
k 2 dkdVL). From (112113ft and (O it follows that the trans- 
fer functions R\ 1+2 (k) do not depend on time variable t. 
That results from the stationarity of the gravitational wave 
field. 

The expressions for the transfer function are calculated 
in the similar way as in (jBaskaran et al . 2008), but now they 
are slightly more complicated: 



k 2 r 2 r 2 
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k 4 c 4 
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fc 2 c 4 cos (kD) sin (kD)kc 4 
4DP + AIP ' 



(21) 



(22) 



These transfer functions behave like k 4 and k 6 respec- 
tively when k — > 0. 

The statistical properties of stochastic gravitational 
wave field may b e characterized by the density parameter 
Q.g W \ Alien! fl996T l. 0, gw is related to the power spectrum 
P h (k): 



Qgw (&) 



2n 2 



Ph(k) 



(23) 



where ku = 2nfn /c = 2-wHo / c, and Ho is the current Hub- 
ble parameter. The density parameter Q gw is the current 
day ratio of energy density of gravitational waves (per unit 
logarithmic interval in k) to the critical density of the Uni- 
verse pcrit = 3c 2 Hq /8tvG. For numerical estimations, we set 
Hubble parameter Ho = 75 — /Mpc and assume a simple 
power law spectrum for the density parameter Q gm : 

( k \ ™ T 

n gw (k) = Q. gw (k ) (j-) ■ (24) 



This form of spectrum can be used as a good approx- 
imation for a large variety of models in frequency range of 
our interest. The flat, scale invariant power spectrum (also 
known as Harrison-Zeldovich power spectrum) corresponds 
to tit = 0. To obtain mean square deviations, we should 
integrate (JT7J) and (O with GW spectrum [f2"3|l. The limits 
of integration fe m i n and fc max are determined from following 
considerations: the highest frequency that puts in the effect 
is defined by cut-off scale coming from the pulsar timing 
technique and determined by total time span of pulsar obser- 
vations used to obtain rotational parameters, the lowest fre- 
quency of GW that can be probed with that method comes 
from the limitations of pulsar-Earth distanc^E &min ~ 75- 
and fc max = 27r/cT b s (note that our maximal frequency 
coincides with minimal frequency of usual pulsar timing 
searches for GWs; in fact fc max = 2an / cT ^ s , a ~ 1, for 
the sake of simplicity we used a — 1). After substitutions 
we obtain the following equations: 
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\ v 
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(26) 



4ir J (n T +2) (^°)^H^0 T [^max 

^n gw (k )k 2 k 2 H ln(>^), n T =2 



where two last terms in eqns. (|21I22|I are omitted be- 
cause they are always much smaller than the first one. More 
stringent constraints come from expression (|26[l because the 
GW influence on the first derivative is hidden inside spin- 
down part. 

As observed value of the second derivative i> bs can not 
exceed the value produced by the effect we can substitute 

< (^) 2 > with ) 2 in ([26|PI . Finally, we arrive at the 
following constraint (adopting the flat spectrum of GW): 



a, 



8ty 2 



V Vo 



T 



2n 2 H 2 \ vo 
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(27) 



(28) 



1 The method can probe cosmological GW with length up to 
present-day horizon size, but effectively their contribution is sup- 
pressed by transfer functions that approach quickly; they should 
only be taken into account in case of GWB spectrum that reddens 
towards lower frequencies. 

2 As a matter of fact, second derivative in pulsar timing analy- 
sis results from fitting a cubic polynomial to observational data. 
Obtained value slightly differs from actual second derivative of 
pulsar residuals owing to the contributions to fitting procedure 
from higher than third order terms and noise terms; I am grate- 
ful to referee for pointing that out. 
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It is instructive to make some numerical estimates us- 
ing the data for PSR B1937+21; that pulsar has been timed 
for a long time with a very high pr ecision. The values of ro - 
tational parameters is taken from (|Manchester et al1l2005t ): 



v = 641 Hz, = 6.2 



n gw < 2.8 x 10" 



Tobs 

10 yrs 



(29) 



these constrai nts are 2 orders of magnitude stronger than the 
current limit (jKopeikinll 19971 ) . It is intriguing that virtually 
all pulsars with small second derivatives exhibit approxi- 
mately the same value of defining relation 10~ 29 < | 1 < 
10~ 28 s -2 0. That relation holds both for ordinary and mil- 
lisecond pulsars. The part in v that is induced by the ef- 
fect, has a T~^ a dependance on total observation time span. 
Having long high-quality timing series we would be able to 
obtain function i'obs(Tobs) and then extract from it the part 
that can be caused by GW influence , thus effectively making 
limits on Q gw tighter. 



3 POSSIBLE METHOD OF DETECTION OF 
ULTRA-LOW FREQUENCY GW 

The technique proposed in the previous section can only 
place limits on the existence of ultra- low frequency GWB; it 
can not be used to detect the presence of that background. 
However, having time series of some parameter that is af- 
fected by the GWB for several pulsars, we can use pairwise 
cross-correla tions in the same way they are used in re sidu- 
als analysis l|Hellings fc Downslll983l ; |jenet et al.ll2005l ). Ac- 
tually, the time se ries for v(t) can be used instead of r(t) 
in "classical" case (|Kaspi. Taylor fc Rvballl994l ; IJenet et al.l 
2006). Data are prepared as follows: total span of timing 
observations is sampled into shorter sub-intervals, e.g. one 
year long and frequency is calculated for every interval. Mag- 
nitude of frequency will decrease linearly due to the regu- 
lar spin-down effect. The first-order contribution from the 
GWB would be absorbed by that much more substantial ef- 
fect. However, the GWB induced effect can be sought in the 
series of observed frequencies after removing of linear trend. 
The correlation coefficient between the observed values of 
Av for each pair of pulsar has the form: 



/(*) 



i= N— 1 

1 Ai/i (t h ei) At/jj(ti,e 2 ) 



V 53 

i=0 
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(30) 



where N - number of sub-intervals, ex, &i - unit vector in 
the directions to pulsars, voi - average values of pulsars ro- 
tational frequencies. 

Due to the linearity, differentiation with respect to 
time does not affect the angular dependence of correla- 
tion function on the angular distance between two pulsars 
(cos (9 = ei ■ e 2 ) and: 



< f(o) >= «tLcw, 



Except PSR J1952+1410 (-£ 



(31) 
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10 



-30 . 



I.6-10" 



- 2 ), J1946+1805 
5.1 ■ 10- 30 s- 2 ). 



- 2 ) and J1823+0550 (-£■ 
Placing that values into J28j produces even more stringent limits 
on U ow . 



c(fl) = 3(1 - cosg) log l- cost 



1 — cos# 1 l.,„, 
~ T" + 2 + 2 



4 2 

where SO equals 1 when 6 — and otherwise. Correlated 
part Av allows us to detect presence of GWB, a\ v =< 

(-) 2 >- 

Correlated part Av induced by GWB can be very 
roughly estimated as: 

Av _ 1 8u(t) (T st x 2 

VQ 



2 v 



2 J 



(32) 



where T si 



length of sub-interval, i.e. one year; on the 



other hand, we can detect fractional deviations of frequency 
4r at R = 10" 14 level. 

Assuming flat spectrum of GWB (Q gw (k) = const) and 
using (|26p we arrive at: 



_Tb_ 

Tsami 



(33) 



where Th = -g- is Hubble time; detection limit is Q, gm ~ 
10~ r . This estimation is certainly too superficial and will be 
reconsidered with use of real pulsar array data. 



4 CONCLUSIONS 

Measuring rotational parameters of pulsar can provide valu- 
able information about processes with characteristic times 
larger than time span of pulsar observations. One of the 
most interesting kind of these processes is influence of ultra 
long-wavelength GWB on pulsar timing. The second deriva- 
tive of rotational frequency can be used to put limits on the 
density of GWB in the frequency region 10~ 12 to 10" 8 Hz. 
The limits for Q gw that come from the suggested method: 

n gw < 2 x io~ 6 . 

Also, the time series of the frequency v from different pulsars 
in pulsar array can be cross-correlated pairwise. Comparison 
of angular dependence of obtained functions with known 
sample of GWB induced correlation similarly to usual pulsar 
timing residual analysis makes possible to detect GWB in 
the frequency range of our interest. That method can probe 
GWB down to energy density Q. gw < 10~ 7 . 
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